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Abstract 

In this paper, we demonstrate that when the ratio n of the number of antenna elements N to the 
number P of radiating sources is superior or equal to 2, then it is possible to choose a propagator from 
a set of n{n + l)/2 — 1 operators to compute the Angles of Arrival (AoA) of the narrowband incoming 
waves. 

This new non eigenbased approach is efficient when the Signal to Noise Ratio (SNR) is moderate, and 
gives multitude of possibilities, that are dependent of the random data, to construct the complex sets 
whose columns are orthogonal to the signal subspace generated by the radiating sources. Elementary 
examples are given for n = 3, n = 4 and n = 6. The simulation results are presented to illustrate the 
performance of the proposed computational methods. 

Keywords: Propagator, DoA, Narrowband, large array, angular spectrum, subspace projector. High 
resolution. 


1 Introduction 

The new technology of planar arrays and nanoarrays consists of constructing thousands of radiating elements 
connected to single unit for processing the emitted and the received signals. With this new conception, the 
computational complexity [1] is proportional to the increasing rate of the antenna elements. When dealing 
with the problem of Angles of Arrival (AoA) detection [1], computing the second order statistics is sufficient, 
however the eigen-approach is computationally heavy. Recent studies focused on the asymptotic behavior 
of the eigenvalues of the cross correlation matrix when its dimensions tend to infinity [2], the distribution 
of the eigenvalues gives the repartition that can be clustered to obtain complex signal and noise subspaces, 
with the latter being characterized by spherical covariance matrix [1]. 

Besides the eigen-approach, other methods [3] were derived by exploiting the evolution of the wavefield 
along the antenna sensors, a linear operator relates the channel matrix of the first antenna elements to the 
rest of them, this operator is called Propagator [4] and can be extracted from the inter spectral matrix. Since 
the appearence of the propagator, few versions were proposed to generalize the concept, in [5] the channel 
matrix was descritized into several partitions where the partial propagator has non redundent elements while 
all the sensors are exploited, in [6], the propagator was extended into third order when the number of sensors 
is superior to two times the number of sources. 

In this paper, we propose a generalization of the propagator for large antenna arrays, we demonstrate 
that given the number of sensor N and the number of emitting sources P, if the ratio n = E{N/P) is 
superior or equal to 2, then it is possible to generate n operators of the noise subspace and the sum of the 
computational methods is given by n{n -I- l)/2 — 1 possibilities. In the next section we present the statistical 
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data model of the problem, in third section we develop the extended propagator and in the last section we 
present some numerical results. 


2 STATISTICAL DATA MODEL 


In this section we elaborate the physical description of the model used in this paper. An array antenna, 
consisting of N identical and isotropic sensors, is receiving a wavefield that is generated by P sources such 
that theirs frequency spectra are concentrated around a carrier frequency fc- We assume that the velocities 
of the sources are negligible during the acquisition time T of the array, dxi/dt ~ 0 for i = 1,..., P, the sources 
are located in the far field region of the antenna. The medium of propagation has the same permettivity 
and permeability (e, fx) in all directions which permits the wavefronts of the sources to be summed linearly. 
If we consider that the noise in the medium, the electronic and thermal noises of the sensors are modeled by 
ergodic random processes, the received signals at instant tk can be written in the following form : 

N P 

X{tk) = EE aj{0i)s^{tk)+ n{tk) (1) 

i=i i=l 


with k = 1,2 ,..., K and K being the number of samples, a{9i) G is the steering vector of the angle 

of arrival, the progressive phase along the sensors is given by : 


a(0i) 


1 p-jhi 

i , C- , C. , . . . , C. 


( 2 ) 


(.)^ is the transposition operator and fXi = 27rdX ^ sinOi, d is the inter-element distance of the array which 
we consider to be uniform and equal to the half of the wavelength A. Si{t) is the carried signal of the 
source with s{t) = [si(t), S 2 (t), sp(t)] and n{t) = [ni{t),n 2 {t), ...,nN{t)] is the additive noise modeled 
by stationary and ergodic zero mean complex random process. In compact form, the received signals are 
given by : 

X{t) = A{9)s{t) + n{t) (3) 

A{9) is the channel matrix of the array. Based on the matrix X{t) we can extract a lot of characteristics of 
the studied system, such as the waveforms of the signals s{t), the polarization of the sources, the coherence 
between the signals s{t) and the number of the signals P. 

In this paper we are focused on detecting the Angles of Arrival (AoA) of the punctual emitters {9^}, for this 
purpose, many methods can be applied to estimate the angles, almost all the techniques are based on second 
order statistics of X {t) which can be computed by : 


r = 


lim 

K—¥-\-0C 


1 

K 


K 


Y,Xit)X+{t) 


(4) 


With finite number of samples K only an approximation is obtained. The theoretical expression of F is given 
by the following equation : 

r = A(0)r,A+(0) + r„ (5) 

the matrix F^ G is the covariance of waveforms, F^ = a^lN is the covariance of the random vector 

n{t), (.)■*■ is the conjugate transpose operator, is the noise variance considered to have the same value 
for all sensors and /at is the identity matrix of dimension N. The subspace based techniques [7-8] exploit 
the orthogonality between the vectorial spaces generated by the signal and noise eigenvectors, indeed the 
spectral theorem is applied to the covariance matrix, so that it can be decomposed into the following form : 

N 

r = Y, = Us<l>sUt + Un^uU+ ( 6 ) 
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Us G is the vectorial space whose columns span the signal subspace and [/„ G £^nxN-p jg gg|. 

whose columns span the noise subspace. The diagonal matrices are given by = diag{Xi,Xp} 

which contains the largest eigenvalues in decreasing order Ai > A 2 > ... > Xp, = diag{Xp+i,X n} 

where the diagonal elements are the noise eigenvalues, in the noiseless case we have Ap+i = ... = Xn = 0. 
The two orthogonal subspaces are related by the following equation : 

UsU+ + UnU+ = In (7) 


In compact form, the spectral decomposition can be written as : 


r= Us Un 




[Us Ur," 


( 8 ) 


With the spectral decomposition, the degenerescence d of the eigenvalues is determined in order to detect 
the number of sources P = N — d hy which the subspace C/„ is delimited to compute the angular spectrum 

[7] : 

/(0,) = {a+(0,)C/„C/+a(0,)} (9) 

for testing angle 9i. If the number P is known, the subspace can be computed using the propagator 
function [4], this operator is described using the channel matrix, theoretically A{6) can be divided into two 
blocks : 

All 


A = 


A2 


( 10 ) 


Ai G and A 2 G given that the two blocks are linearly dependent, a linear operator 11 G 

(j^N-PxP gxist such that A 2 = AAi and 11 = A 2 A^^, consequently the vectorial space whose columns span 
the noise subspace C/„ is given by 

[n I —In-p] [A{9)] = On-pxp (11) 

This operator is extracted from the covariance matrix, similarly P is split into two blocks T = [Tg | r„] with 
Tg G , r„ G an estimate of the propagator is given by the following equation : 

n+ = rir„ (12) 

with (.)! is the Penrose pseudo inverse, r| = (r+rs)“^r+. If we consider the noiseless case, using the 
equation (10), the covariance matrix is decomposed into four blocks : 

Til ri2 

r2i I'22 


T = 


(13) 


with Til G r 2 i G ri 2 G and r 22 G ^ Following the decomposition of 

the channel matrix [5], the F blocks are computed by : r 22 = A 2 TsA", Fh = AiF^yl^ and F 21 = A 2 TsA". 
We can obeserve that F 21 = IIFu which leads to a noise subspace : 

[r2iFj“j^^ I —In-p] [A{9)] = On-pxp (14) 

Usually in the litterature [4],[9-14] the propagator is developped using the equation (12), the second 
approach in equation (14) was proposed in [5], however we can deduce another version based on the equation 
nl = (A 2 A'^^)" the corresponding subspace is described by the following equation : 


-Ip I ri2Fi2 


[A(9)\ = Opxp 


(15) 


All the three versions are equivalent, but in the presence of the matrix n(t) and the variation of the numbers 
N and P, they are different due to the presence of the noise powers in the diagonal blocks Fn and r 22 . 
Given this situation, if we are supposed to work with a new architecture of the arrays comprising of hundreds 
of radiating elements which might be the case for nanoarrays, the question that we can ask is given a large 
number N ( N » P) can we find other possibilities for computing the propagator operator without using 
the diagonal blocks Fj^ ? the answer of this question is dependent on the ratio n = N/P which will give a 
choice to choose an operator from a set that we will elaborate in the next section. 
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3 Extended Propagator 


In this section, we develop our approach based on Uniform Linear Array (ULA) which is considered to be 
centro-symmetric, taking the first element of the array as the origin of coordinates, the channel matrix has 
a Vandermonde structure. 


A = 


( 




1 \ 

e-JMP 




g-i(iV-l)MP j 


(16) 


where A has the properties Aq A* = Inxp, rank{A) = P and JmA* = AAp where Jn is the exchange 
matrix given by : 


/ 

0 

0 ... 

0 

1 \ 


0 

0 ... 

1 

0 


: 0 

1 ... 

: 0 

: 0 

V 

1 

0 ... 

0 

oy 


Ap is P X P diagonal matrix Ap = (iiag{e we note that the rank of A is dependent 

on the azimuths of the sources 9i. Let us consider the three approaches for the propagator presented earlier: 


Q-- 

Qi 


-- '(rtr„)+ I -i^.p 
= [r2ir;^]^^ I —In-p_ 


Q 2 — 


-Ip I ri2rL 


(17) 


The maximum number of sources for Q to resolve is V — 1, which is also the case for Qi, however if V < 2P, 
Q 2 will fail to compute the null space of A{9), this problem is due to the density of matrix Q 2 which 
contains P(P — 1) zeros, if we compute the value of spectrum with arbitrary vector a{9i) then the function 
a{9i)Q2Q2aA{9i) will not give us information whether 9i is an Angle Of Incidence (AOI) or not. 

Supposing that the Signal to Noise Ratio (SNR) is high, evaluating Q 2 requires the condition N/P > 2, 
following the same reasoning the ratio must also satisfy n < [V/PJ, so the variable n is bounded : 

2 < n < L^J (18) 

We introduce the selection matrices [16] to make the formalism rigourous, let us consider the set {ei}i<i<„ 
where a matrix Ci G is defined by Cj = [Opxp(i-i) | Ip \ OpxN-{Pi)\'^ which verifies the following 

equations : 

e[ej=S^jIp (19) 


Y^e,eJ = lN ( 20 ) 

n 

'^eje^di = Ip ( 21 ) 

\\e,\\F = VP (22) 

where di = 1/n and P denotes the Frobenius norm (||e||p = y^Tr{ee~^)). The channel matrix is decomposed 
into : 

n n 

A = A = (23) 

2=1 2=1 
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This decomposition is equivalent to ^ the blocks Ai of rank F are equivalent of scalars 

in Euclidean base with the following property: 

A = efA, A+ = A+e, (24) 

This partition is mapped into partitions on the covariance matrix T, this concept is the same principle 
used in the semiclassical approach of statistical mechanics where the one dimenisonal space phase is divided 
into cells on the condition that the product of the step position and impulsion equals the quantum constant 
h, in our case the division of T implies dXdY = this concept is illustrated in Fig.l. 
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Figure 1: Similarity concept of partition principle between statistical physics and array processing : (a) 
semiclassical approach of one dimensional space phase (p, q), (b) uniform partition of matrix T where the 
mark + denotes the affected blocks by the noise Tu and N = 12P. 


Given n phases of A there exist n? — n linear operators Qij that relates the and phases Ai = 
QijAj. The reverse time evolution of the wavefronts is introduced {i < j) unlike the standard propagator 
{i = 2 > j = 1). The microstate between the and the sets of the sensors is given by : 

Tij = ef Tej = ef ATsA+e^ = AiTsA^ (25) 


fii — + cAlp 


(26) 


Given the spectral matrix with the property T’*' = T the corresponding conjugate transpose block T)^- is 
computed by : 


r 


+ 


(efrej)~^ = eJVci = Tji 


(27) 


For infinite number of samples K, the general partition of T is : 


/ fn 

r2i 


ri2 

f22 


Tm \ 
r2n 


V 







Where r„„ 
by: 


AnTgA^ + from equation (25), the partitions of the channel matrix are related 


efA = TifeT^^ejA 
Ai = TifeT^^Tl 
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This equation proves that we have (n — 2) possibilities for arbitrary index k ^ i ^ j, k can be chosen 
randomly because we assume the existence of isotropic noise r„ = u^In- Using equation (28) we can derive 
the following equation : 

n 

r,fcr]^ejA-(n-l)efA = OpxP (29) 

The above equation allows us to compute the N x P matrix whose columns span the null space of the channel 
matrix A as follows : 

\ E r,fcr],ej -/3pef U = OpxP (30) 

Where (3p = {n — l)/p, the above equation is the propagator based on the block Ai, hence we have 
n possibilities which can be concatenated in single matrix G as follows : 


-Pp 

Qi2 

Qi3 

Qln 

Q 21 

-/3p 

Q 23 ■■■ 

Q2n 

Qnl 

Qn2 

Qn3 

—/?7V-(n-l)P / 


(31) 


Given two blocks (i, j) the propagator that maps the phase of the wavefront from block i to the corresponding 
j is product of partial propagators of blocks i < k < j in both causal case i > j and retarded case i < j: 


1-1 

Qij — Qk^k — 1 

k—i 


i-1 

Qij — Qk,k-\-l 

k—i 


(32) 


Ik has the property Tr(iE') = —(n — 1)N and each block iky = ikjj : 


= (eJ’^'Ci)^ (33) 

To illustrate the general structure of ik, let us take an example of uniform linear array consisting of TV = 3 
antennas and one radiating source P = 1, the number of possibilities of partition is n = {2, 3}, let us take 
n = 3, the orthonormal base is reduced into {ei}i<i <3 where [ei, 62 , 63 ] = I3. The steering vector is given by 

A = lei+e^‘^ie 2 + e2i‘^ie3 (34) 


with (j)i = 2 'k\ sin 01 , applying the previous formalism yields to the following result of ik : 


1 

to 


e-2jy.i 

4' = 

1 

to 



eJd>i 

-2 


Hence the rows of ik are orthogonal to the steering vector A : 

^^A = ^2A = 4'3H = 0 


(35) 


(36) 


We can deduce that iky = ik^/ for single source P = 1, ik is self adjoint operator ik+ = ik. The general rule 
elaborated is that for every integer n, there are partitions on T and n possible propagators are derived. 
Given sufficiently high SNR, we can state the following proposition : 

Proposition 3.1. Given a random matrix X G £^xk ^ ^ CA/’(0Arxi, T), 3n G N such that V2 < n < 
, 3 H G i = 1, ..., n} with Card{^) = n{n + l)/2 — 1 where we have : 


lim ik„H = Oatxp 

SNR ^+00 


(37) 
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(38) 


IfP=l : 


with the conditions k ^ j ^ i and 13p 


lim '^niA = Opxp 
SNR^+oo 

n — 1 

'Hm = Y, r,fcr]^ej - /3p 
i=i 


eT 


'^N = 


(n - l)/p 


(39) 

(40) 


If we take an example of an array consisting of iV = 500 cells and five statistically independent sources 
with moderate SNR and sufficient number of snapshots we have card{fl) = *) “ 5049 possibilities 

n = {'I' 2 ) 4 ^ 3 ,'I'loo} where the limit partition is n = 100. To make the idea more clear, we present in the 
next section an example for n = 3. 


4 Elementary example, n = 3 

In this example, we suppose that n = 3 = [N/P\, the constants become f3p = 2Jp, i = {1,2,3} and the 
generated base is {ei}i<i< 3 , the channel matrix is decomposed into the following : 

A = cic^ A 62^2 A 6363 A (41) 


— siAi + 62212 + 63213 

The dimensions of the three blocks are 2 I 1 G A 2 G and 2 I 3 G C^- 22 ^xp g^g . 



(42) 


The partial operators are given by : 

2I1 = Ql2^2j2li = Qi3A3,A2 = Q23A3,Qi3 = Q 12 Q 23 
Consequently the spectral matrix is partitioned into the following : 

r = 


'll 

ri2 

ri3 

'21 

f22 

r23 

31 

r32 

f33 


(43) 


(44) 


The index k has n — 2 = 1 possibility, let us compute the first propagator : 

vl/3i = - 2 ef+ Y 

j=2,k=^j 


— —26^ + ri3r23e^ + ri2r32e3 

Similarly we obtain the two left propagators as follows : 

4'32 = r23r|3e{’ — 2e^ + r2ir3;^e3 
4'33 = r32r|2cf + r3ir2ie^ — 263 


(45) 
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(46) 

(47) 



The assembled matrix ^'3 is given by : 


^3 


-21 p 


ri242 \ 


r23rl3 

-2/p 

r2irL 

(48) 

r32rl2 

r3irL 

—2In-2P / 



A particular case 'I'aa in this example was presented in [ 6 ], where it was stated that the method is comparable 
to the standard Propagator and standard MUSIC estimators when the SNR is sufficiently high, besides it 
was also stated that ^'33 can resolve two closely spatial sources if the noise is spatially non uniform. The 
other operators ^'31 and 4*32 are not mentioned in the literature. 


5 summary 

In this brieve section, we describe the algorithm presented earlier : 


Table I: Extended Propagator algorithm. 


Input : (r,P). 

1. Choose n from 2 < n < IN/P\. 

2 . Generate {ei}i<i<„ G 

3. Choose p from {1, 2,..., n}. 

d.Compute rpfcPjfceJ - Ppe^-jk ^ j ^ p). 


In the following table we present different cases encountred in real life applications and the cases where 
the proposed formalism is available : 


0 < n < 1 

D = 0 , possibility exists for quasi-stationary signals using KHATRI RAO space [15]. 

1 < n < 2 

D = 0, Standard Propagator is available. 

2 < n < IN/P\ 

^ = {'I' 2 , n{n - 1 - l )/2 — 1 possibilities of extended propagators. 


The proposed operators are valid for arbitrary geometry as long as the distance between consecutive 
sensors is about half the wavelength : 

||u+i-r'||<^ (49) 

With being the position of the sensor per given reference. 

The accuracy of the spectra is measured by comparing theirs variance with the Cramer Rao Bound (CRB). 
But in terms of localization function, we propose the following definition of the ideal spectrum. The one 
dimensional ideal spectrum f{6) for detecting P sources is Dirac comb : 


lim 

SNR^+00 


p 

f{e) = Y,ak5{e-e^) 


(50) 


with 0 G other formulation, the ideal spectrum is the one sided Fourier transform of superpo¬ 

sition of P sinusoidals sampled with frequency F > tt : 


p 

k^l 
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(51) 












With u!k = ^irOk and q = 9~^, most of the subspace based techniques do not directly provide the information 
about the magnitudes ak, the obtained peaks measure the degree of orthogonality between the signal and 
noise subspaces. In the next section we present some simulation results to test the performance of several 
propagators 'i’np comparatively to some high resolution eigen-based methods such as ESPRIT which is briefly 
described in the following section. 

5.1 Basic Review of ESPRIT algorithm 

The ESPRIT algorithm ( Estimation of Signal Parameters via Rotational Invariance) [8]-[16] is high resolu¬ 
tion and double eigendecomposition approach, it is based on computing the DoAs from the eigenvalues of 
rotation operator that relates two identical partitions of the antenna array given that each element in one 
subarray has its own twin in the other, in other words, the array is characterized by displacement invariance 
d. The follwoing table summarizes the algorithm : 

Table 2: ESPRIT algorithm. 


Input : X 

1. Create two partitions of array by choosing the number m defined by: P < m < N — 1. 

2. Compute signal subspace Us from X = [/EC+ or from P = UAU~^. 

3. Generate selection matrices ei = [im \ Omx(Ar-m)] € e 2 = [Omx(Af-m) I Im] and solve : 

eiUs'S ~ e2Us. 

4. Compute the eigenvalues of I* = T^T~^. 

5. Extract the DoAs by = Arg{^i} and 9i = Arcsin{—Xfj,i/2TTd) for I < i < P 


6 SIMULATION RESULTS 

We run some computer simulations to verify the performance of some propagators presented in this paper, 
for this purpose we consider P = 3 narrowband punctual sources with carrier frequency of fc = 2GHz 
are impinging on an array consisting of = 18 omnidirectional sensors uniformely spaced with half the 
wavelength d = 7cm which makes the total length of the array Lx = 119cm with corresponding Rayleigh 
angular limit resolution of 9hpbw — 6 °. 

The sources are located in the Fraunhofer region where theirs azimuths, measured from broadside, are 
9i = 10°, 02 = 21° and 03 = 45°. The received signals are equipowered and complex identically distributed 
random processes Si{t) ~ CAf{0,<Ts,i = Iwatt) and the number of snapshots is A" = 200 samples . 

The number of possibilities is given by the ensemble = { 4 > 2 , ^' 3 , 'I' 4 , di'g} with Card{Lt) = 20, due to 
the diversity of possibilities we only take some operators for performance evaluation. 

In the first experiment, we fix SNR = 5dB, the figure 2 represents an average of P = 200 Monte Carlo runs 
of the operators, 'I' 2 i,'I' 3 i,'I '32 and 'I' 4 i. 

The peaks are sharp for all the operators and the physical interpretation is easy to make due to lack of 
sidelobs, however the difference lies in the magnitude of the peaks for each 4'^. 

Using the same conditions, we compute the correlation coefficients of 4*4 to explore the relationship between 
four operators, the following matrix represents the coefficients : 


/ 1.00 

0.87 

0.80 

0.97 \ 

0.87 

1.00 

0.96 

0.87 

0.80 

0.96 

1.00 

0.79 

0.97 

0.87 

0.79 

1.00 / 


We observe that the pairs ( 4 > 4 i, 'I' 44 ) and ( 4 ' 43 , 'I' 42 ) are highly correlated. The figure 3 represents the Root 
Mean Square Error (RMSE) of 4*4 comparatively to the standard ESPRIT method [ 8 ]. When SNR is about 
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Figure 2: Average oi L = 200 Monte Carlo runs of {'£' 21 , 'I'ai, '£' 32 , '£' 41 } witlr d = A/2, SNR = 5dB, K = 200 
and 0 = [ 10°, 21°, 45° ]. 



Figure 3: Average of L = 200 Monte Carlo runs of 4^4 and ESPRIT for each value of SNR, d = A/2, 
K = 200 and 0 = [ 10°, 21°, 45° ]. 


20dB, all the spectra are equivalent, in the other hand we realize that ('£' 41 ,'£' 44 ) vary with the same rate, 
which is also the same case for ( 4 ' 43 , '£' 42 ), this remark is confirmed by the matrix 

£n the last experiment, we compare the spectrum of with that MUSIC [7] and G-MUS£C [2], the last 
method outperforms the two other while the function of d'gi presents several fluctuations as minor sidelobs. 

Some of the immediate enhancements that can be added to the present paper is the orthonormalization 
[4] of any propagator '^np which will add more resolution power in source detection. The other proposition 
lies in using the real valued transformation [17] of the spectral matrix F, this will reduce the computational 
complexity and increase the probability of locating the coherent sources. 

In the other hand, the drawback of the proposed formalism resides in noise condition, the SNR must be 
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Figure 4: Average oi L = 200 Monte Carlo runs of three spectra with SNR = lOdB, d = A/2, K = 200 and 
e = [ 10°, 21°, 45° ]. 


moderate and the number of radiating sources P must be known, hence a non eigen-based spectral method 
to determine P will complete the extended propagator algorithm. 

7 Conclusion 

In this paper, we have introduced a proposition which states that when the ratio n of the number of antenna 
sensors N to the number of narrowband non coherent sources P is superior or equal to 2, and the SNR 
is sufficiently high, it is possible to generate n{n -I- l)/2 — 1 distinct versions of high resolution propagator 
operators to detect the angular positions of the radiating sources with no constraint on the geometry of the 
array. 

Original formalism is presented, inspired from semiclassical approach of statistical mechanics, to give the limit 
of the channel matrix partition, from which diversity of extended propagators is derived, this new approach 
is suitable for large array. In Monte Carlo simulation results, we have presented elementary examples to test 
the performance of the extended propagators of orders n = 3, n = 4 and n = 6 which proved theirs resolving 
powers in high SNR conditions. 


8 Appendix 

In this section we present parts of computer programs implemented in simulations, using MATLAB software 

• Operators 'I'ai and 4'32 : 

N=I8; 

P=3; 

RI2=R(I:P,P-fI:N); 

RI3=R(I:P,2*P-fI:N); 

R23=R(P-bI:2*P,2*P-bI:N); 

RI2=R(I:P,P-fI:2*P); 

R2I=RI2’; 

R32=R23’; 
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R31=R13’; 

Q31=[-2*eye(P) R13*pinv(R23) R12*pinv(R32)]; 
Q32=[R23*pinv(R13) -2*eye(P) R21*pinv(R31)]; 


• Operator '£'41 : 

R13=R(1:P,2*P+1:3*P); 

R23=R(P+1:2*P,2*P+1:3*P); 

R34=R(2*P+1:3*P,3*P+1:N); 

R14=R(1:P,3*P+1:N); 

R42=R(3*P+1:N,P+1:2*P); 

R12=R(1:P,P+1:2*P); 

Q41=[-3*eye(P) R13*pinv(R23) R14*pinv(R34) R12*pmv(R42)]; 


• Operator 'I'ei : 

R12=R(1:P,P+1:2*P); 

R13=R(1:P,2*P+1:3*P); 

R23=R(P+1:2*P,2*P+1:3*P); 

R14=R(1:P,3*P+1:4*P); 

R34=R(2*P+1:3*P,3*P+1:4*P); 

R15=R(1:P,4*P+1:5*P); 

R45=R(3*P+1:4*P,4*P+1:5*P); 

R16=R(1:P,5*P+1:N); 

R56=R(4*P+1:5*P,5*P+1:N); 

R63=R(5*P+1:N,2*P+1:3*P); 

Q61=[-5*eye(P) R13*pinv(R23) R14*pinv(R34) R15*pinv(R45) R16*pinv(R56) R13*pinv(R63)]; 
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